Dependencies and data

# dependencies
library(tidyverse)
library(knitr)
library(kableExtra)
library(boot)
library(parallel)
library(bayestestR)
library(patchwork)
library(mdthemes)
library(lme4)
library(sjPlot)
library(emmeans)
library(ggstance)
library(janitor)
# library(merTools) called via merTools:: to avoid namespace collisions between MASS and dplyr


# set seed for reproducibility
set.seed(42)

# options
options(knitr.table.format = "html") # necessary configuration of tables

# disable scientific notation
options(scipen = 999) 

# function to round all numeric vars in a data frame
round_df <- function(df, n_digits = 3) {
  require(janitor)
  df %>% mutate_if(is.numeric, janitor::round_half_up, digits = n_digits)
}

# create necessary directories
dir.create("../data/processed")
dir.create("../data/results")
#dir.create("models")

# get data 
# data with confidence intervals
data_estimates_D <- read_csv("../data/processed/data_estimates_D.csv") %>%
  filter(method == "bca")

data_estimates_iat_D <- read_csv("../data/processed/data_estimates_iat_D.csv") %>%
  mutate(trial_type = "iat",
         unique_id = as.factor(unique_id))

data_demographics_iat <- read_rds("../data/processed/data_iat_processed_participant_level.rds") %>%
  mutate(session_id = as.factor(session_id)) %>%
  semi_join(data_estimates_iat_D, by = c("session_id" = "unique_id")) %>%
  select(unique_id = session_id, age, sex)

Sample descriptives

data_demographics_iat %>%
  summarize(min_age  = round_half_up(min(age, na.rm = TRUE), 2),
            max_age  = round_half_up(max(age, na.rm = TRUE), 2),
            mean_age = round_half_up(mean(age, na.rm = TRUE), 2),
            sd_age   = round_half_up(sd(age, na.rm = TRUE), 2)) %>%
  gather() %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
key value
min_age 8.00
max_age 81.00
mean_age 31.60
sd_age 12.46 
data_demographics_iat %>%
  count(sex) %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
sex n
f 6296
m 3193
NA 11

CI widths

MAP-MAP width

Most probable estimate among the most probable estimates

data_map_ci_widths <- data_estimates_iat_D %>%
  group_by(domain, trial_type) %>%
  do(point_estimate(.$ci_width, centrality = "MAP")) %>%
  ungroup()

write_csv(data_map_ci_widths, "../data/results/data_map_ci_widths_iat_d.csv") 

data_map_ci_widths %>%
  summarize(map_map = point_estimate(MAP, centrality = "MAP"),
            min_map = min(MAP),
            max_map = max(MAP)) %>%
  unnest(map_map) %>%
  rename(MAP_MAP = MAP) %>%
  round_df(2) %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
MAP_MAP min_map max_map
0.71 0.47 0.72

MAP CI width

By domain and trial type

data_map_ci_widths %>%
  pivot_wider(names_from = trial_type, values_from = MAP) %>%
  round_df(2) %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
domain iat
50 Cent - Britney Spears 0.70
African Americans - European Americans 0.70
Artists - Musicians 0.69
Asians - Whites 0.71
Astrology - Science 0.70
Atheism - Religion 0.71
Athletic People - Intelligent People 0.72
Avoiding - Approaching 0.51
Bill Clinton - Hillary Clinton 0.70
Briefs - Boxers 0.69
Burger King - McDonald’s 0.70
Canadian - American 0.71
Capital Punishment - Imprisonment 0.71
Career - Family 0.70
Chaos - Order 0.56
Coffee - Tea 0.70
Cold - Hot 0.71
Conservatives - Liberals 0.68
Corporations - Nonprofits 0.71
David Letterman - Jay Leno 0.71
Denzel Washington - Tom Cruise 0.69
Determinism - Free will 0.67
Difficult - Simple 0.66
Dogs - Cats 0.66
Dramas - Comedies 0.69
Drinking - Abstaining 0.71
Effort - Talent 0.67
Evolution - Creationism 0.70
Fat People - Thin People 0.70
Foreign Places - American Places 0.71
Friends - Family 0.71
Gay People - Straight People 0.71
George Bush - John Kerry 0.71
Gun Control - Gun Rights 0.71
Helpers - Leaders 0.70
Hiphop - Classical 0.71
Innocence - Wisdom 0.70
Japan - United States 0.71
Jazz - Teen Pop 0.70
Jews - Christians 0.71
Jocks - Nerds 0.72
Kobe - Shaq 0.71
Lawyers - Politicians 0.71
Lord of the Rings - Harry Potter 0.70
Manufactured - Natural 0.63
Meat - Vegetables 0.69
Meg Ryan - Julia Roberts 0.70
Microsoft - Apple 0.71
Money - Love 0.47
Mother Teresa - Princess Diana 0.71
Mountains - Ocean 0.71
Muslims - Jews 0.70
National Defense - Education 0.69
New York - California 0.71
Night - Morning 0.68
Numbers - Letters 0.70
Old People - Young People 0.71
Organized Labor - Management 0.71
Pants - Skirts 0.71
Past - Future 0.65
Pepsi - Coke 0.71
Poor People - Rich People 0.56
Private - Public 0.70
Prolife - Prochoice 0.71
Protein - Carbohydrates 0.71
Protestants - Catholics 0.71
Punishment - Forgiveness 0.57
Realism - Idealism 0.71
Reason - Emotions 0.71
Rebellious - Conforming 0.68
Receiving - Giving 0.71
Redsox - Yankees 0.71
Relaxing - Exercising 0.71
Republicans - Democrats 0.68
Rich People - Beautiful People 0.71
Security - Freedom 0.71
Single - Married 0.71
Skeptical - Trusting 0.51
Solitude - Companionship 0.70
Southerners - Northerners 0.70
Speed - Accuracy 0.71
Stable - Flexible 0.71
State - Church 0.71
Strong - Sensitive 0.71
Tall People - Short People 0.70
Tax Reductions - Social Programs 0.71
Team - Individual 0.71
Technology - Nature 0.69
Television - Books 0.70
Tradition - Progress 0.71
Traditional Values - Feminism 0.71
Urban - Rural 0.71
West Coast - East Coast 0.71
Winter - Summer 0.71
Wrinkles - Plastic Surgery 0.71

Plot by domain and trial type

data_ci_width_map_D <- data_estimates_iat_D %>%
  group_by(domain, trial_type) %>%
  do(point_estimate(.$ci_width, centrality = "MAP")) %>%
  ungroup() %>%
  mutate(MAP = round_half_up(MAP, 3),
         trial_type = case_when(trial_type == "tt1" ~ "Trial type 1",
                                trial_type == "tt2" ~ "Trial type 2",
                                trial_type == "tt3" ~ "Trial type 3",
                                trial_type == "tt4" ~ "Trial type 4",
                                trial_type == "iat" ~ "IAT"),
         trial_type = fct_relevel(trial_type, "Trial type 1", "Trial type 2", "Trial type 3", "Trial type 4", "IAT")) %>%
  mutate(domain = fct_rev(domain))

# save to disk
write_csv(data_ci_width_map_D, "../data/results/data_ci_width_map_iat_D.csv")

# plot
p_ci_widths <- 
  ggplot(data_ci_width_map_D, aes(MAP, domain)) + 
  geom_point(position = position_dodge(width = 0.8)) +
  mdthemes::md_theme_linedraw() +
  #facet_wrap(~ trial_type, ncol = 4, nrow = 1) +
  labs(x = "Highest probability (MAP) 95% CI width",
       y = "") + 
  theme(legend.position = "top")

p_ci_widths

Proportion different from zero

Caterpillar plot

p_cis_by_domain <- 
  data_estimates_iat_D %>%
  arrange(estimate) %>%
  group_by(domain) %>%
  mutate(ordered_id = row_number()/n()) %>%
  ungroup() %>%
  ggplot() +
  geom_linerange(aes(x = ordered_id, ymin = ci_lower, ymax = ci_upper, color = sig),
                 alpha = 1) +
  geom_point(aes(ordered_id, estimate), size = 0.5, shape = "square") +
  geom_hline(yintercept = 0, linetype = "dotted") +
  mdthemes::md_theme_linedraw() +
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        legend.position = "top") +
  scale_color_viridis_d(end = 0.6, direction = -1) +
  xlab("Ranked participant") +
  ylab("*D* score") +
  labs(color = "95% CI excludes zero point") + 
  facet_wrap(~ domain, ncol = 10)

p_cis_by_domain

Calculate scores

data_diff_zero <- 
  bind_rows(
    mutate(data_estimates_D, measure = "IRAP"),
    mutate(data_estimates_iat_D, measure = "IAT")
  ) %>%
  mutate(measure = fct_relevel(measure, "IRAP", "IAT"),
         domain = as.factor(domain),
         trial_type = case_when(trial_type == "tt1" ~ "Trial type 1",
                                trial_type == "tt2" ~ "Trial type 2",
                                trial_type == "tt3" ~ "Trial type 3",
                                trial_type == "tt4" ~ "Trial type 4",
                                trial_type == "iat" ~ "IAT"),
         trial_type = fct_relevel(trial_type, "Trial type 1", "Trial type 2", "Trial type 3", "Trial type 4", "IAT")) %>%
  group_by(domain, trial_type, measure) %>%
  summarize(proportion_diff_zero = mean(sig),
            variance = plotrix::std.error(sig)^2,
            .groups = "drop") %>%
  # model cannot be run on zero variance or 0 or 1 logit, so offset by a minuscule amount
  mutate(proportion_diff_zero_temp = case_when(proportion_diff_zero < 0.001 ~ 0.001, 
                                               proportion_diff_zero > 0.999 ~ 0.999,
                                               TRUE ~ proportion_diff_zero),
         proportion_diff_zero_logit = boot::logit(proportion_diff_zero_temp)) %>%
  select(-proportion_diff_zero_temp) %>%
  #filter(!(proportion_diff_zero == 0 & variance == 0)) %>%
  mutate(variance = ifelse(variance == 0, 0.0001, variance)) 

# data_diff_zero %>%
#   round_df(2) %>%
#   kable() %>%
#   kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)

# save to disk
write_csv(data_diff_zero, "../data/results/data_diff_zero_irap_vs_iat.csv")

Plot

p_diff_zero <- 
  data_diff_zero %>%
  filter(measure == "IAT") %>%
  mutate(domain = fct_rev(factor(domain))) %>%
  ggplot(aes(proportion_diff_zero, domain)) +
  geom_linerangeh(aes(xmin = proportion_diff_zero - sqrt(variance)*1.96,
                      xmax = proportion_diff_zero + sqrt(variance)*1.96),
                  position = position_dodge(width = 0.75)) + 
  geom_point(position = position_dodge(width = 0.75)) +
  #scale_shape_manual(labels = c("IRAP", "IAT"), values = c(15, 16)) +
  #scale_color_viridis_d(begin = 0.3, end = 0.7, labels = c("IRAP", "IAT")) +
  mdthemes::md_theme_linedraw() +
  labs(x = "Proportion of scores<br/>different from zero point",
       y = "") + 
  theme(legend.position = "top",
        panel.spacing = unit(1.5, "lines")) +
  coord_cartesian(xlim = c(0,1))

p_diff_zero

Model

NB model is slightly different to the one used to compare IRAP D and PI scores: that one has (1) no random slope for measure and (2) a random intercept for trial type too. Including (1) seemed important given that the two IRAP and IAT demonstrated very different heterogeneity between domains. Not including it greatly and inappropriately expands the prediction intervals on the IRAP (i.e., variation in the IAT is modeled as variation in both, in appropriately). In contrast, the effects were very similar between IRAP D and IRAP PI, so this wasn’t necessary in the other analysis. Including (2) gave convergence issues, likely because the IAT only has a single trial type, so it was dropped.

# fit model
fit_diff_zero <- 
  lmer(proportion_diff_zero_logit ~ 1 + measure + (measure | domain),
       weights = 1/variance, 
       data = data_diff_zero,
       # solution from https://www.metafor-project.org/doku.php/tips:rma_vs_lm_lme_lmer
       control = lmerControl(check.nobs.vs.nlev = "ignore",  
                             check.nobs.vs.nRE = "ignore"))

# extract marginal means
results_emm_diff_zero <- 
  summary(emmeans(fit_diff_zero, ~ measure)) %>%
  dplyr::select(measure, estimate = emmean, se = SE, ci_lower = lower.CL, ci_upper = upper.CL)

# extract re Tau
results_re_tau_diff_zero <- fit_diff_zero %>%
  merTools::REsdExtract() %>%
  as_tibble(rownames = "re") %>%
  rename(tau = value) 

# combine
results_diff_zero <- results_emm_diff_zero %>%
  mutate(pi_lower = estimate - (1.96 * sqrt(se^2 + results_re_tau_diff_zero$tau^2)),  # as in metafor package's implementation of prediction intervals, see metafor::predict.rma.R 
         pi_upper = estimate + (1.96 * sqrt(se^2 + results_re_tau_diff_zero$tau^2))) |>
  select(-se) |>
  mutate_if(is.numeric, boot::inv.logit)

# plot
p_prop_nonzero <- 
  ggplot(results_diff_zero, aes(measure, estimate)) +
  geom_linerange(aes(ymin = pi_lower, ymax = pi_upper), size = 0.5, position = position_dodge(width = 0.8), linetype = "dotted") +
  geom_linerange(aes(ymin = ci_lower, ymax = ci_upper), size = 1.25, position = position_dodge(width = 0.8)) +
  geom_point(position = position_dodge(width = 0.8), size = 2.5) +
  mdthemes::md_theme_linedraw() +
  scale_y_continuous(breaks = c(0, .25, .5, .75, 1), labels = c("0.00<br/>(Worse)", "0.25", "0.50", "0.75", "1.00<br/>(Better)")) +
  #scale_color_viridis_d(alpha = 1, begin = 0.3, end = 0.7, labels = c("IRAP", "IAT")) +
  #scale_shape_manual(labels = c("IRAP", "IAT"), values = c(15, 16)) +
  scale_x_discrete(labels = c("IRAP D scores", "IAT D scores")) +
  labs(x = "",
       y = "Proportion of scores<br/>different from zero point<br/>") + 
  theme(legend.position = "none") +
  coord_flip(ylim = c(0, 1))

p_prop_nonzero

results_diff_zero %>%
  round_df(2) %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
measure estimate ci_lower ci_upper pi_lower pi_upper
IRAP 0.06 0.04 0.10 0.00 0.62
IAT 0.57 0.52 0.62 0.06 0.97 
# tests
data_emms_diff_zero <- emmeans(fit_diff_zero, list(pairwise ~ measure), adjust = "holm") 

summary(data_emms_diff_zero)$`pairwise differences of measure` %>%
  as.data.frame() %>%
  select(comparison = 1, p.value) %>%
  mutate(p.value = ifelse(p.value < .001, "< .001", round_half_up(p.value, 3))) %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
comparison p.value
IRAP - IAT < .001

Proportion different from one another

Within domain and trial type.

Note: Discriminability between a score and zero can be determined using the CI, because zero is a known value and only the score is measured with uncertainty. However, discriminability between two scores must take into account the uncertainty in the estimation of both scores. Weir (2005) argues that such an interval can be estimated by expanding the CIs by sqrt(2). Here I refer to these intervals as Discriminability Intervals (DIs).

Calculate discriminability

# # discriminability using non-overlap of CIs
# discriminability <- function(data, i) {
#   data_with_indexes <- data[i,] # boot function requires data and index
#   ci_lower <- data_with_indexes$ci_lower 
#   ci_upper <- data_with_indexes$ci_upper
#   n_ci_lower <- length(ci_lower)
#   n_ci_upper <- length(ci_upper)
#   r_ci_lower <- sum(rank(c(ci_lower, ci_upper))[1:n_ci_lower])
#   A <- (r_ci_lower / n_ci_lower - (n_ci_lower + 1) / 2) / n_ci_upper
#   return(A)
# }

# discriminatory using the significance of the difference score
# the goal here is to assess mean_diff > 1.96 * sqrt(SE1^2 + SE2^2 for every possible comparison EXCLUDING self comparisons. This is tricky to do within a typical tidyverse workflow as it means doing mutates involving each row of a column and every other row of that column but not the same row.
  # the below solution is to use expand.grid to find all combinations of a row with itself, and then use the modulus of the length of the row to filter out the self-pairings. Then do mutates on the rows to assess significant differences. It's enough to then summarize the proportion of significant results across all participants.
discriminability <- function(data, i) {
  data_with_indexes <- data[i,] # boot function requires data and index
  
  grid_estimates <- expand.grid(data_with_indexes$estimate, data_with_indexes$estimate) |>
    mutate(diff = Var1 - Var2,
           row_number = row_number(),
           modulus = row_number %% (nrow(data_with_indexes)+1)) |>
    filter(modulus != 1) |>
    select(diff)
  
  grid_se <- expand.grid(data_with_indexes$se, data_with_indexes$se) |>
    mutate(critical_value = 1.96 * sqrt(Var1^2 + Var2^2),
           row_number = row_number(),
           modulus = row_number %% (nrow(data_with_indexes)+1)) |>
    filter(modulus != 1) |>
    select(critical_value)
  
  proportion_sig_diff <- 
    bind_cols(grid_estimates, grid_se) |>
    mutate(sig = abs(diff) > critical_value) |>
    summarize(proportion_sig_diff = mean(sig)) |>
    pull(proportion_sig_diff)
  
  return(proportion_sig_diff)
}

bootstrap_discriminability <- function(data){
  
  require(dplyr)
  require(boot)
  
  fit <- 
    boot::boot(data      = data, 
               statistic = discriminability, 
               R         = 2000,
               sim       = "ordinary", 
               stype     = "i",
               parallel  = "multicore", 
               ncpus     = parallel::detectCores()-1)
  
  results <- boot::boot.ci(fit, conf = 0.95, type = "bca") 
  
  output <-
    tibble(
      estimate = fit$t0,
      ci_lower = results$bca[4],
      ci_upper = results$bca[5]
    )
  
  return(output)
}

# irap data
data_discriminability_D <- read_csv("../data/results/data_discriminability_D.csv") %>%
    filter(method == "bca")

# bootstrapping has a long execution time, so load saved values if they've already been calculated
if(file.exists("../data/results/data_discriminability_iat_D.csv")) {
  
  data_discriminability_iat_D <- read_csv("../data/results/data_discriminability_iat_D.csv")
  
} else {
  
  # bootstrap D scores 
  data_discriminability_iat_D <- data_estimates_iat_D |>
    mutate(se = (ci_upper - ci_lower)/(1.96*2)) |>
    select(unique_id, domain, trial_type, estimate, se) |>
    group_by(domain, trial_type) |>
    do(bootstrap_discriminability(data = .)) |>
    ungroup() |>
    rename(proportion_discriminable = estimate) |>
    mutate(variance = (((ci_upper - ci_lower)/(1.96*2)))^2,
           domain = as.factor(domain),
           #trial_type = fct_relevel(trial_type, "tt1", "tt2", "tt3", "tt4", "iat"),
           measure = "IAT") 
  
  # save to disk
  write_csv(data_discriminability_iat_D, "../data/results/data_discriminability_iat_D.csv")
  
}

Plot

# combine
data_discriminability_combined <- 
  bind_rows(
    mutate(data_discriminability_D, measure = "IRAP"),
    mutate(data_discriminability_iat_D, measure = "IAT")
  ) %>%
  mutate(measure = fct_relevel(measure, "IRAP", "IAT"),
         trial_type = case_when(trial_type == "tt1" ~ "Trial type 1",
                                trial_type == "tt2" ~ "Trial type 2",   
                                trial_type == "tt3" ~ "Trial type 3",   
                                trial_type == "tt4" ~ "Trial type 4",
                                trial_type == "iat" ~ "IAT"),
         trial_type = fct_relevel(trial_type, "Trial type 1", "Trial type 2", "Trial type 3", "Trial type 4", "IAT")) %>%
  #filter(!(proportion_discriminable == 0 & variance == 0)) %>%
  mutate(variance = ifelse(variance == 0, 0.0001, variance)) |>
  # model cannot be run on zero variance or 0 or 1 logit, so offset by a minuscule amount
  mutate(
    proportion_discriminable_temp = case_when(proportion_discriminable < 0.001 ~ 0.001, 
                                              proportion_discriminable > 0.999 ~ 0.999,
                                              TRUE ~ proportion_discriminable),
    proportion_discriminable_logit = boot::logit(proportion_discriminable_temp)
  ) %>%
  select(-proportion_discriminable_temp)

p_discriminability <- 
  data_discriminability_combined %>%
  filter(measure == "IAT") %>%
  mutate(domain = fct_rev(factor(domain))) %>%
  ggplot(aes(proportion_discriminable, domain)) +
  geom_linerangeh(aes(xmin = proportion_discriminable - sqrt(variance)*1.96,
                      xmax = proportion_discriminable + sqrt(variance)*1.96),
                  position = position_dodge(width = 0.75)) + 
  geom_point(position = position_dodge(width = 0.75)) +
  #scale_shape_manual(labels = c("IRAP", "IAT"), values = c(15, 16)) +
  #scale_color_viridis_d(begin = 0.3, end = 0.7, labels = c("IRAP", "IAT")) +
  mdthemes::md_theme_linedraw() +
  #facet_wrap(~ trial_type, ncol = 4) +
  labs(x =  "Proportion of scores<br/>differerent from other scores",
       y = "") + 
  theme(legend.position = "top",
        panel.spacing = unit(1.5, "lines")) +
  coord_cartesian(xlim = c(0,1))

p_discriminability

Model

# fit meta analytic model
fit_disciminability <- 
  lmer(proportion_discriminable_logit ~ 1 + measure + (measure | domain), 
       weights = 1/variance, 
       data = data_discriminability_combined,
       # solution from https://www.metafor-project.org/doku.php/tips:rma_vs_lm_lme_lmer
       control = lmerControl(check.nobs.vs.nlev = "ignore",  
                             check.nobs.vs.nRE = "ignore"))

# extract marginal means
results_emm_disciminability <-
  summary(emmeans(fit_disciminability, ~ measure)) %>%
  dplyr::select(measure, estimate = emmean, se = SE, ci_lower = lower.CL, ci_upper = upper.CL) 

# extract re Tau
results_re_tau_disciminability <- fit_disciminability %>%
  merTools::REsdExtract() %>%
  as_tibble(rownames = "re") %>%
  rename(tau = value) 

# combine
results_disciminability <- results_emm_disciminability %>%
  mutate(pi_lower = estimate - (1.96 * sqrt(se^2 + results_re_tau_disciminability$tau^2)),  # as in metafor package's implementation of credibility intervals, see metafor::predict.rma.R
         pi_upper = estimate + (1.96 * sqrt(se^2 + results_re_tau_disciminability$tau^2))) |>
  select(-se) |>
  mutate_if(is.numeric, boot::inv.logit)

# plot
p_prop_discriminable <-
  ggplot(results_disciminability, aes(measure, estimate)) +
  geom_linerange(aes(ymin = pi_lower, ymax = pi_upper), size = 0.5, position = position_dodge(width = 0.8), linetype = "dotted") +
  geom_linerange(aes(ymin = ci_lower, ymax = ci_upper), size = 1.25, position = position_dodge(width = 0.8)) +
  geom_point(position = position_dodge(width = 0.8), size = 2.5) +
  scale_y_continuous(breaks = c(0, .25, .5, .75, 1), labels = c("0.00<br/>(Worse)", "0.25", "0.50", "0.75", "1.00<br/>(Better)")) +
  #scale_shape_manual(labels = c("IRAP", "IAT"), values = c(15, 16)) +
  #scale_color_viridis_d(begin = 0.3, end = 0.7, labels = c("IRAP", "IAT")) +
  scale_x_discrete(labels = c("IRAP D scores", "IAT D scores")) +
  mdthemes::md_theme_linedraw() +
  labs(x = "",
       y = "Proportion of scores<br/>differerent from other scores<br/>") +
  theme(legend.position = "none") +
  coord_flip(ylim = c(0, 1))

p_prop_discriminable 

results_disciminability %>%
  round_df(2) %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
measure estimate ci_lower ci_upper pi_lower pi_upper
IRAP 0.07 0.05 0.10 0.01 0.47
IAT 0.46 0.42 0.49 0.07 0.91 
# tests
data_emms_disciminability <- emmeans(fit_disciminability, list(pairwise ~ measure), adjust = "holm")

summary(data_emms_disciminability)$`pairwise differences of measure` %>%
  as.data.frame() %>%
  select(comparison = 1, p.value) %>%
  mutate(p.value = ifelse(p.value < .001, "< .001", round_half_up(p.value, 3))) %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
comparison p.value
IRAP - IAT < .001

CI widths as a proportion of observed range

NB observed range of confidence intervals

Calculate scores

## calculate observed ranges 
observed_range_estimates_D <- data_estimates_D %>%
  group_by(domain, trial_type) %>%
  dplyr::summarize(min = min(ci_lower, na.rm = TRUE),
                   max = max(ci_upper, na.rm = TRUE),
                   .groups = "drop") %>%
  mutate(range = max - min) 

observed_range_estimates_iat_D <- data_estimates_iat_D %>%
  group_by(domain) %>%
  dplyr::summarize(min = min(ci_lower, na.rm = TRUE),
                   max = max(ci_upper, na.rm = TRUE),
                   .groups = "drop") %>%
  mutate(range = max - min) 

# calculate CI / range 
data_ci_width_proportions_D <- data_estimates_D %>%
  # join this data into the original data
  full_join(observed_range_estimates_D, by = c("domain", "trial_type")) %>%
  # calculate ci width as a proportion of observed range
  mutate(ci_width_proportion = ci_width / range) %>%
  mutate(measure = "IRAP") 

data_ci_width_proportions_iat_D <- data_estimates_iat_D %>%
  # join this data into the original data
  full_join(observed_range_estimates_iat_D, by = "domain") %>%
  # calculate ci width as a proportion of observed range
  mutate(ci_width_proportion = ci_width / range) %>%
  mutate(measure = "IAT")

# combine
data_ci_width_proportions_combined <- 
  bind_rows(
    data_ci_width_proportions_D,
    data_ci_width_proportions_iat_D
  ) %>%
  mutate(measure = fct_relevel(measure, "IRAP", "IAT"),
         domain = as.factor(domain),
         trial_type = fct_relevel(trial_type, "tt1", "tt2", "tt3", "tt4", "iat")) %>%
  group_by(measure, domain, trial_type) %>%
  summarize(ci_width_proportion_mean = mean(ci_width_proportion),
            variance = plotrix::std.error(ci_width_proportion)^2) %>%
  ungroup() %>%
  # logit transform
  mutate(ci_width_proportion_mean_temp = case_when(ci_width_proportion_mean < 0.0001 ~ 0.0001,
                                                   ci_width_proportion_mean > 0.9999 ~ 0.9999,
                                                   TRUE ~ ci_width_proportion_mean),
         ci_width_proportion_mean_logit = boot::logit(ci_width_proportion_mean_temp)) %>%
  select(-ci_width_proportion_mean_temp)

write_csv(data_ci_width_proportions_combined, "../data/results/data_ci_width_proportions_irap_d_vs_iat_d.csv")

Plot

p_coverage <- 
  data_ci_width_proportions_combined %>%
  mutate(domain = fct_rev(factor(domain))) %>%
  filter(measure == "IAT") %>%
  ggplot(aes(ci_width_proportion_mean, domain)) +
  geom_point(position = position_dodge(width = 0.75)) +
  scale_shape_manual(labels = c("*D* scores", "PI scores"), values = c(15, 16)) +
  geom_linerangeh(aes(xmin = ci_width_proportion_mean - sqrt(variance)*1.96,
                      xmax = ci_width_proportion_mean + sqrt(variance)*1.96),
                  position = position_dodge(width = 0.75)) + 
  scale_color_viridis_d(begin = 0.3, end = 0.7, labels = c("*D* scores", "PI scores")) +
  mdthemes::md_theme_linedraw() +
  #facet_wrap(~ trial_type, ncol = 4) +
  labs(x = "Proportion of observed range covered<br/>by individual scores' 95% CIs",
       y = "") + 
  theme(legend.position = "top",
        panel.spacing = unit(1.5, "lines")) +
  coord_cartesian(xlim = c(0,1))

p_coverage

Model

# fit model
fit_ci_width_proportions <- 
  lmer(ci_width_proportion_mean_logit ~ 1 + measure + (measure | domain), 
       weights = 1/variance,
       data = data_ci_width_proportions_combined,
       # solution from https://www.metafor-project.org/doku.php/tips:rma_vs_lm_lme_lmer
       control = lmerControl(check.nobs.vs.nlev = "ignore",  
                             check.nobs.vs.nRE = "ignore"))

# extract marginal means
results_emm_ci_width_proportions <-
  summary(emmeans(fit_ci_width_proportions, ~ measure)) %>%
  dplyr::select(measure, estimate = emmean, se = SE, ci_lower = lower.CL, ci_upper = upper.CL)

# extract re Tau
results_re_tau_ci_width_proportions <-
  merTools::REsdExtract(fit_ci_width_proportions) %>%
  as_tibble(rownames = "re") %>%
  rename(tau = value)

# combine
results_ci_width_proportions <- results_emm_ci_width_proportions %>%
  mutate(pi_lower = estimate - (1.96 * sqrt(se^2 + results_re_tau_ci_width_proportions$tau^2)),  # as in metafor package's implementation of credibility intervals, see metafor::predict.rma.R
         pi_upper = estimate + (1.96 * sqrt(se^2 + results_re_tau_ci_width_proportions$tau^2))) %>%
  select(-se) %>%
  mutate_if(is.numeric, boot::inv.logit)

# plot
p_ci_width_proportion_observed_range <-
  ggplot(results_ci_width_proportions, aes(measure, estimate, 
  )) +
  geom_linerange(aes(ymin = pi_lower, ymax = pi_upper), size = 0.5, position = position_dodge(width = 0.8), linetype = "dotted") +
  geom_linerange(aes(ymin = ci_lower, ymax = ci_upper), size = 1.25, position = position_dodge(width = 0.8)) +
  geom_point(position = position_dodge(width = 0.8), size = 2.5) +
  #scale_shape_discrete(labels = c("IRAP", "IAT")) +
  scale_y_continuous(breaks = c(0, .25, .5, .75, 1), labels = c("0.00<br/>(Better)", "0.25", "0.50", "0.75", "1.00<br/>(Worse)")) +
  #scale_shape_manual(labels = c("IRAP", "IAT"), values = c(15, 16)) +
  #scale_color_viridis_d(begin = 0.3, end = 0.7, labels = c("IRAP", "IAT")) +
  scale_x_discrete(labels = c("IRAP D scores", "IAT D scores")) +
  mdthemes::md_theme_linedraw() +
  labs(x = "",
       y = "Proportion of observed range covered<br/>by individual scores' 95% CIs") +
  theme(legend.position = "none") +
  coord_flip(ylim = c(0, 1))

p_ci_width_proportion_observed_range

results_ci_width_proportions %>%
  round_df(2) %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
measure estimate ci_lower ci_upper pi_lower pi_upper
IRAP 0.51 0.49 0.53 0.42 0.60
IAT 0.27 0.26 0.27 0.16 0.41 
# tests
data_emms_ci_width_proportions <- emmeans(fit_ci_width_proportions, list(pairwise ~ measure), adjust = "holm") 

summary(data_emms_ci_width_proportions)$`pairwise differences of measure` %>%
  as.data.frame() %>%
  select(comparison = 1, p.value) %>%
  mutate(p.value = ifelse(p.value < .001, "< .001", round_half_up(p.value, 3))) %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
comparison p.value
IRAP - IAT < .001

Combined plots

Supplementary Figure 2S

p_cis_by_domain

ggsave(filename  = "plots/supplementary_figure_2S_cis_by_domain_iat_d.pdf",
       plot      = p_cis_by_domain,
       device    = "pdf",
       # path      = NULL,
       # dpi       = 300,
       units     = "in",
       width     = 16,
       height    = 16,
       limitsize = TRUE)

Supplementary Figure 3S

p_ci_widths

ggsave(filename  = "plots/supplementary_figure_3S_ci_widths_iat_d.pdf",
       plot      = p_ci_widths,
       device    = "pdf",
       # path      = NULL,
       # dpi       = 300,
       units     = "in",
       width     = 6,
       height    = 12,
       limitsize = TRUE)

Supplementary Figure 4S

p_diff_zero

ggsave(filename  = "plots/supplementary_figure_4S_proportion_excluding_zero_point_iat.pdf",
       plot      = p_diff_zero,
       device    = "pdf",
       # path      = NULL,
       # dpi       = 300,
       units     = "in",
       width     = 6,
       height    = 14,
       limitsize = TRUE)

Supplementary Figure 5S

p_discriminability

ggsave(filename  = "plots/supplementary_figure_5S_proportion_discriminable_iat.pdf",
       plot      = p_discriminability,
       device    = "pdf",
       # path      = NULL,
       # dpi       = 300,
       units     = "in",
       width     = 6,
       height    = 14,
       limitsize = TRUE)

Supplementary Figure 6S

p_coverage

ggsave(filename  = "plots/supplementary_figure_6S_proportion_coverage_iat.pdf",
       plot      = p_coverage,
       device    = "pdf",
       # path      = NULL,
       # dpi       = 300,
       units     = "in",
       width     = 6,
       height    = 14,
       limitsize = TRUE)

Figure 6

p_combined <- 
  p_prop_nonzero + 
  p_prop_discriminable + 
  p_ci_width_proportion_observed_range +
  plot_layout(ncol = 1) #, guides = "collect") & theme(legend.position = "top")

p_combined

ggsave(filename  = "plots/figure_6_metaanalyses_irap_d_vs_iat_d.pdf",
       plot      = p_combined,
       device    = "pdf",
       # path      = NULL,
       # dpi       = 300,
       units     = "in",
       width     = 5,
       height    = 5,
       limitsize = TRUE)

Session info

sessionInfo()
## R version 4.2.1 (2022-06-23)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur ... 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] janitor_2.1.0     ggstance_0.3.5    emmeans_1.7.5     sjPlot_2.8.10    
##  [5] lme4_1.1-30       Matrix_1.4-1      mdthemes_0.1.0    patchwork_1.1.1  
##  [9] bayestestR_0.12.1 boot_1.3-28       kableExtra_1.3.4  knitr_1.39       
## [13] forcats_0.5.1     stringr_1.4.0     dplyr_1.0.9       purrr_0.3.4      
## [17] readr_2.1.2       tidyr_1.2.0       tibble_3.1.8      ggplot2_3.3.6    
## [21] tidyverse_1.3.2  
## 
## loaded via a namespace (and not attached):
##   [1] TH.data_1.1-1       googledrive_2.0.0   minqa_1.2.4        
##   [4] colorspace_2.0-3    ellipsis_0.3.2      sjlabelled_1.2.0   
##   [7] estimability_1.4    snakecase_0.11.0    markdown_1.1       
##  [10] parameters_0.18.1   fs_1.5.2            gridtext_0.1.4     
##  [13] ggtext_0.1.1        rstudioapi_0.13     listenv_0.8.0      
##  [16] furrr_0.3.0         farver_2.1.1        bit64_4.0.5        
##  [19] fansi_1.0.3         mvtnorm_1.1-3       lubridate_1.8.0    
##  [22] xml2_1.3.3          codetools_0.2-18    splines_4.2.1      
##  [25] cachem_1.0.6        sjmisc_2.8.9        jsonlite_1.8.0     
##  [28] nloptr_2.0.3        ggeffects_1.1.2     pbkrtest_0.5.1     
##  [31] broom_1.0.0         dbplyr_2.2.1        broom.mixed_0.2.9.4
##  [34] shiny_1.7.2         effectsize_0.7.0    compiler_4.2.1     
##  [37] httr_1.4.3          sjstats_0.18.1      backports_1.4.1    
##  [40] assertthat_0.2.1    fastmap_1.1.0       gargle_1.2.0       
##  [43] cli_3.3.0           later_1.3.0         htmltools_0.5.3    
##  [46] tools_4.2.1         coda_0.19-4         gtable_0.3.0       
##  [49] glue_1.6.2          merTools_0.5.2      Rcpp_1.0.9         
##  [52] cellranger_1.1.0    jquerylib_0.1.4     vctrs_0.4.1        
##  [55] svglite_2.1.0       nlme_3.1-157        iterators_1.0.14   
##  [58] insight_0.18.0      xfun_0.31           globals_0.15.1     
##  [61] rvest_1.0.2         mime_0.12           lifecycle_1.0.1    
##  [64] googlesheets4_1.0.0 future_1.27.0       MASS_7.3-57        
##  [67] zoo_1.8-10          scales_1.2.0        vroom_1.5.7        
##  [70] promises_1.2.0.1    hms_1.1.1           sandwich_3.0-2     
##  [73] yaml_2.3.5          sass_0.4.2          stringi_1.7.8      
##  [76] highr_0.9           foreach_1.5.2       plotrix_3.8-2      
##  [79] blme_1.0-5          rlang_1.0.4         pkgconfig_2.0.3    
##  [82] systemfonts_1.0.4   arm_1.12-2          evaluate_0.15      
##  [85] lattice_0.20-45     labeling_0.4.2      bit_4.0.4          
##  [88] tidyselect_1.1.2    parallelly_1.32.1   magrittr_2.0.3     
##  [91] R6_2.5.1            generics_0.1.3      multcomp_1.4-20    
##  [94] DBI_1.1.3           pillar_1.8.0        haven_2.5.0        
##  [97] withr_2.5.0         abind_1.4-5         survival_3.3-1     
## [100] datawizard_0.4.1    performance_0.9.1   modelr_0.1.8       
## [103] crayon_1.5.1        utf8_1.2.2          tzdb_0.3.0         
## [106] rmarkdown_2.14      grid_4.2.1          readxl_1.4.0       
## [109] reprex_2.0.1        digest_0.6.29       webshot_0.5.3      
## [112] xtable_1.8-4        httpuv_1.6.5        munsell_0.5.0      
## [115] viridisLite_0.4.0   bslib_0.4.0